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D ■ Abstract 

^^ , The abundance of collapsed objects in the universe, or halo mass function, is an important theoretical tool 
0\J ' in studying the effects of primordially generated non-Gaussianities on the large scale structure. The non- 
Gaussian mass function has been calculated by several authors in different ways, typically by exploiting the 
smallness of certain parameters which naturally appear in the calculation, to set up a perturbative expansion. 
We improve upon the existing results for the mass function by combining path integral methods and saddle 
point techniques (which have been separately applied in previous approaches). Additionally, we carefully 
i-^ ' account for the various scale dependent combinations of small parameters which appear. Some of these com- 
binations in fact become of order unity for large mass scales and at high redshifts, and must therefore be 
O ' treated non-perturbatively. Our approach allows us to do this, and also to account for multi-scale density 
4-J , correlations which appear in the calculation. We thus derive an accurate expression for the mass function 
jrt I which is based on approximations that are valid over a larger range of mass scales and redshifts than those 
of other authors. By tracking the terms ignored in the analysis, we estimate theoretical errors for our result 
and also for the results of others. We also discuss the complications introduced by the choice of smoothing 
^ ' filter function, which we take to be a top-hat in real space, and which leads to the dominant errors in our 
CO , expression. Finally, we present a detailed comparison between the various expressions for the mass functions, 
^r\ ' exploring the accuracy and range of validity of each. 
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O ! 1 Introduction 

The primordial curvature inhomogeneities, generated by the inflationary mechanism, obey nearly Gaussian 
1 J - statistics. The deviations from Gaussianity, while expected to be small, provide a unique window into the 
i'jjj , physics of inflation. For example, single-field slow-roll models of inflation lead to a small level of non- 
C^ ' Gaussianity (NG), so that an observation of a large NG would indicate a deviation from this paradigm. 

Until a few years ago, the main tool to constrain NG was considered to be the statistics of the cosmic 
microwave background (CMB) temperature field, since inhomogeneities at the CMB epoch are small and their 
physics can be described by a perturbative treatment. In recent years, however, thanks to observations and 
developments in the theory, the large-scale structure (LSS) of the universe has emerged as a complementary 
probe to constrain primordial NG. While it is true that the n-point functions of the density field on small scales 
are dominated by the recent gravitational evolution, and do not reflect anymore the statistics of primordial 
perturbations, it turns out that the abundance of very massive objects, which form out of high peaks of the 
density perturbations, is a powerful probe of primordial NG. In this context, much attention has been given 
recently to three possible methods of constraining the magnitude and shape of the primordial NG with the LSS: 
the galaxy power spectrum, the galaxy bispectrum and the mass function. It was pointed out in Refs. [1, 2] 
that a NG of a local type induces a scale dependence on the galaxy power spectrum, thus making it a sensitive 
probe of the magnitude of local NG /j^£. From Ref. [3] one finds the following constraints: —29 < /^£ < -1-69, 
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already comparable with those obtained from CMB measurements in Ref. [4]: — 10 < /^£ < +74. The future 
is even more promising, with precisions of A/^£ ^ 10 [5] and A/^£ ^ 1 [6, 7, 8] being claimed for future 
surveys. The galaxy bispectrum is also a promising probe of NG as it could be more sensitive to other triangle 
configurations [9] . The mass function - which is the focus of this work, and which we discuss in detail below 
- has been used for example in Ref. [5] together with the scale dependent bias to produce forecasts for future 
surveys, and in Ref. [10] in an attempt to explain the presence of a very massive cluster at a large redshift as 
an indication of a large NG. For more references and information we refer the reader to reviews summarizing 
recent results on these topics [11, 12]. 

The formation of bound dark matter halos from initially small density perturbations, as seen in numerical 
simulations, is a complicated and violent process. Some insight into the physics involved has been gained 
from the study of analytical models. The quantity of interest is the halo mass function, defined as the number 
density of dark matter halos with a mass between M and M + dM, 
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where p is the average density of the universe, a{M) is the variance of the density contrast Sr filtered on 
some comoving scale R corresponding to the mass M , and the function /(ct) is to be computed. Throughout 
this work, we will refer to /(cr) itself as the mass function. A very useful tool in the analysis is the spherical 
collapse model [13], which predicts that the value of the linearly extrapolated density contrast of a spherical 
halo, at the time when the halo collapses, is Sc — 1.686, with a weak cosmology dependence. This value serves 
as a collapse threshold for determining which inhomogeneous regions will end up as collapsed objects. Using 
this idea, Press & Schechter [14] (PS) first computed the mass function /(a) in the case of Gaussian initial 
conditions. Their calculation however suffered from a problem of undcrcounting which affects the overall 
normalization - their approach does not count underdense regions embedded in larger overdense regions 
as eventually collapsed objects. To account for this discrepancy, PS introduced an ad-hoc factor of 2 by 
demanding that the mass function be correctly normalized, such that all the mass in the universe must be 
contained in collapsed objects. In the excursion set approach. Bond et al. [15] resolved this issue and derived 
a correctly normalized mass function, for Gaussian initial conditions. They argued that the filtered density 
contrast 5ji follows a random walk as a function of the filtering scale, and the problem of computing ,f{(j) is 
translated into the problem of finding the rate of "first crossing" of the barrier 5c, whose solution is well-known. 
We will study this formalism in detail in section 3 for the more general non-Gaussian case. 

Turning to non-Gaussianities, the most popular non-Gaussian mass functions are those due to Matarrese, 
Verde and Jimenez [16] (MVJ) and LoVerde et al. [17] (LMSV). Both groups used the PS approach, by 
modifying the probability density function for the (linearized) density contrast to describe non-Gaussian 
initial conditions. In their prescription, the relevant object is the ratio i?ng of non-Gaussian to Gaussian mass 
functions. The full mass function is usually taken as the product of i?ng and an appropriate Gaussian mass 
function as given by iV-body simulations, e.g. the Sheth & Tormen mass function [18]. It is not clear however 
that this is the correct way to proceed. Indeed, in a series of papers [19, 20, 21], Maggiore & Riotto (MR) 
presented a rigorous approach to the first-passage problem in terms of path integrals, and in Ref. [21] they 
pointed out that a PS-like prescription in fact misses some important non-Gaussian effects stemming from 
3-point correlations between different scales (so-called "unequal time" correlators). 

On the other hand, MR treated non-Gaussian contributions to /(cr) by simply linearizing in the 3-point 
function of Sr, i.e. by linearizing in the non-Gaussian parameter /nl- Since the NG are assumed to be 
small, in the sense that the parameter e — {S^)/a^ satisfies e ^ 1, one might expect that such a perturbative 
treatment is valid. However, another crucial ingredient in the problem is that the length scales of interest are 
large, which leads to a second small parameter v"^ where i/ — Sc/cr. This is evident in the calculations of MR, 
who crucially use t^~^ oc cr^ as a small parameter. Any perturbative treatment now depends not only on the 
smallness of e and u^^ individually, but also on the specific combinations of these parameters which appear 
in the calculation. It is known (and we will explicitly see below) that a natural combination that appears 
is eiy^, which can become of order unity on scales of interest. The mass functions given by LMSV and MR 
therefore break down as valid series expansions when this occurs. Interestingly, MVJ's PS-like treatment on 
the other hand involved a saddle point approximation, allowing them to non-perturbatively account for the 
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term (which appears in an exponential in their approach). For a discussion, see Ref. [22]. 



It appears to us therefore, that there is considerable room for improvement in the theoretical calculation 
of the mass function. The goal of our paper is twofold. Firstly, we present a rigorous calculation of the mass 
function in the following way : (a) we use the techniques developed by MR in Refs. [19, 20, 21], which allow us 
to track the complex multi-scale correlations involved in the calculation, and (b) we demonstrate that MR's 
approach can be combined with saddle point techniques (used by MVJ), to non-perturbatively handle terms 
which can become of order unity. This leads to an expression for the mass function which is valid on much 
larger scales than those presented by MR and LMSV. Secondly, by keeping track of the terms ignored, we 
calculate theoretical error bars on the expressions for /(a) resulting not only from our own calculations, but 
also for those of the other authors [16, 17, 21]. Since the terms ignored depend on u in general, these error 
bars are clearly scale dependent. This allows us to estimate the validity of each of the expressions for the 
mass function at different scales, but importantly it also allows us to analytically compare between different 
expressions. In this paper we will not explicitly account for effects of the ellipsoidal collapse model [23, 24], 
since these are expected to be negligible on the very large scales which are of interest to us. For a recent 
treatment of ellipsoidal collapse effects in the presence of non-Gaussianities on scales where ev^ ^ 1, see 
Lam & Sheth [25]. For a different approach to computing the non-Gaussian halo mass function, see Ref. [1], 
where the authors proposed that this mass function can be approximated as a convolution of the Gaussian 
mass function with a probability distribution function that maps between halos identified in Gaussian and 
non-Gaussian A'^-body simulations. This probability distribution itself was approximated as a Gaussian, with 
mean and variance fit from simulations. While in practice this approach is easy to implement, it relies heavily 
on the output of iV-body simulations. Our approach, on the other hand, allows us to compute a mass function 
almost entirely from first principles. 

This paper is organized as follows. In section 2 we fix some notation and briefly introduce the two most 
popular shapes of primordial NG, i.e. the local and equilateral ones. In section 3 we present our calculation 
of the mass function. In section 4 we discuss certain subtleties regarding the truncation of the perturbative 
series, and also compare with the other expressions for /(cr) mentioned above. In section 5 we discuss the 
effects induced by some additional complications introduced in the problem due to the specific choice of the 
filter function [19], which we take to be a top-hat in real space, and due to the inclusion of stochasticity in 
the value of the collapse threshold 5c (which is also expected to partially account for effects of ellipsoidal 
collapse) [20] . In section 6 we compare our final result Eqn. (43) with those of other authors, including 
theoretical errors for each, and conclude with a brief discussion of the results and directions for future work. 
Some technical asides have been relegated to the Appendices. 

2 Models of non-Gaussianity 

We need to relate the linearly evolved density field to the primordial curvature perturbation, which carries the 
information of the non-linearities produced during and after inflation. We start from the Bardeen potential 
$ on subhorizon scales, given by 

$(k,z) = -^T(fc):^7^(fc), (2) 
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where 7^(k) is the (comoving) curvature perturbation, which stays constant on superhorizon scales; T{k) is 
the transfer function of perturbations, normalized to unity as /c — > 0, which describes the suppression of power 
for modes that entered the horizon before the matter-radiation equality; and D{z) is the linear growth factor 
of density fluctuations, normalized such that D{z) = {l-\-z)~^ in the matter dominated era. Then, the density 
contrast field is related to the potential by the Poisson equation, which in Fourier space reads 

<5(k, z) = -,^-^$(k, z) = -^^T{k)D{z)n{k) 

= M{k,z)n{k), (3) 

where we substituted Eqn. (2). Here, i7m is the present time fractional density of matter (cold dark matter 
and baryons), and _ffo = 100/ikms^ Mpc^ is the present time Hubble constant. The redshift dependence 



is trivially accounted for by the linear growth factor D(z) and in the following, for notational simplicity, we 
will often suppress it. All our calculations will use a reference ACDM cosmology compatible with WMAP7 
data [4], using parameters h = 0.702, ftm — 0.272, present baryon density fib = 0.0455, scalar spectral 
index Us = 0.961 and as — 0.809, where erf is the variance of the density field smoothed on a length scale 
of 8/i^^Mpc. For simplicity, for the transfer function T{k) we use the BBKS form, proposed in Bardcen et 
al. [23]: 

rBBKs(a;) = — ^ In (1 + 2.34.T) (l + 3.89a; + (le.la;)^ + {hAQxf + {QJlxfY^'^ , (4) 

l.iAx 

where x = A:(/iMpc~ )/r with a shape parameter T = fim/icxp — 0;,(1 + y2h/il„i) that accounts for bary- 

onic effects as described in Ref. [26]. For more accurate results, one could use a numerical transfer function, 
as obtained by codes like CMBFAST [27] or CAMB [28]; the results are not expected to be qualitatively 
different. 

In order to study halos, which form where an extended region of space has an average overdensity which 
is above threshold, it is useful to introduce a filter function Wij(|x|), and consider the smoothed density field 
(around one point, which we take as the origin), 

/fjSu , 
^-^W{kR)S{k), (5) 

where W{kR) is the Fourier transform of the filter function. For all numerical calculations we will use the 
spherical top-hat filter in real space, whose Fourier transform W{kR) is given by 

3 

W{y) = -^ (sin y-ycosy) . (6) 

This choice allows us to have a well-defined relation between length scales and masses, namely M — (4:Tr/'S)flmPcR 
with pc = 3Hq/{8ttG) — 2.75 • 10^^/i^^Afsoi(^~^Mpc)~'^. However it introduces some complexities in the anal- 
ysis, which we will comment on later. By using Eqns. (5) and (3) we have, for the 3-point function, 

/^3i, j3jL j3i, , , , 

(2^)3(271)3(271)3 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^''^^^^'^^^^^''^^^^' 

. . ^^) 

where the subscript c denotes the connected part, and analogous formulae are valid for the higher order 
correlations. 

2.1 Shapes of non-Gaussianity 

The function (7^(ki)7^(k2)7?.(k3) )c encodes information about the physics of the infiationary epoch. By 
translational invariance, it is proportional to a momentum-conserving delta function: 

(7^(kl)7^(k2)7^(k3) )c - (27r)-''r5i,(ki + k2 + k3)Bn{kuk2, fcg) , (8) 

where the (reduced) bispectrum B-j^^ki, fc2, k^) depends only on the magnitude of the fc's by rotational invari- 
ance. According to the particular model of infiation, the bispectrum will be peaked about a particular shape 
of the triangle. The two most common cases are the squeezed (or local) NG, peaked on squeezed triangles 
ki <^ k2 — ks, and the equilateral NG, peaked on equilateral triangles ki ~ k2 — ks. Indeed, one can define 
a scalar product of bispectra, which describes how sensitive one is to a NG of a given type if the analysis is 
performed using some template form for the bispectrum. As expected, the local and equilateral shapes are 
approximately orthogonal with respect to this scalar product [29] . We will now describe these two models in 
more detail. 

The local model: 
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The local bispcctruni is produced when the NG is generated outside the horizon, for instance in the curvaton 
model [30, 31] or in the inhomogeneous reheating scenario [32]. In these models, the curvature perturbation 
can be written in the following form, 

7^(x) ^ 7e,(x) + |/;°£ (7e^(x) - {UD) + ^.gNL7^^(x) , (9) 

where TZg is the linear, Gaussian field. We have included also a cubic term, which will generate the trispectrum 
at leading order. The bispcctrum is given by 

Bn{ki,k2, h) - ^f^°l [Pn{ki)Pn{k2) + cycl.] , (10) 

5 

where "cycl." denotes the 2 cyclic permutations of the wavenumbers, and Pjiik) is the power spectrum given 
by P-Tzik) = ^fc""'~^. The trispectrum is given by 



;7^(kl)7^(k2)7^(k3)7^(k4) )c - (27r)3,5z,(ki + ka + kg + k4) 
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25/NL E ^^(l''« + k''l)^^(^&)^^(^-)+25-'^NL J2 Pnika)Pnikb)Pnikc) 

b<c a<b<c 



(11) 



The equilateral model: 

Models with derivative interactions of the inflaton field [33, 34, 35] give a bispcctrum which is peaked around 
equilateral configurations, whose specific functional form is model dependent. Moreover, the form of the 
bispcctrum is usually not convenient to use in numerical analyses. This is why, when dealing with equilateral 
NG, it is convenient to use the following parametrization, given in Ref. [36], 
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Bn{kiMM)^^f^i^ A ^,4-„.,4-». + on. L L ^2f4-„^/.s - /,. ,.2,.3^r4-„^/.s +5 perms. . (12) 



2fcf-"=fc^-"= 3(fcifc2A:3)2(4-"e)/3 (fcifc2fc|)(4-"a)/3 



This is peaked on equilateral configurations, and its scalar product with the bispectra produced by the realistic 
models cited above is very close to one. Therefore, being a sum of factorizable functions, it is the standard 
template used in data analyses. 

3 Random walks and the halo mass function 

We now turn to the main calculation of the paper. The non-Gaussian halo mass function can be obtained by 
calculating the barrier first crossing rate J^ of a random walk with non-Gaussian noise, in the presence of an 
absorbing barrier. This can be done perturbatively, starting from a path integral approach as prescribed by 
MR [19, 21] and the mass function can be shown to be /(cr) = 2a'^F{<j). As discussed by MR, the calculation 
of / involves certain assumptions regarding the type of filter used and also the location of the barrier. In 
particular, the formalism is simplest for a sharp filter in fc-space, and using the spherical top-hat of Eqn. 
(6) introduces complications in the form of non-Markovian effects. Further, in order to make the spherical 
collapse ansatz more realistic and obtain better agreement with iV-body simulations, MR show that it is 
useful to treat the location of the barrier 5c as a stochastic variable itself, and allow it to diffuse. For the time 
being, we will ignore these complications, and will return to their effects in section 5. 

To make the paper self-contained, we begin with a brief review of the path integral approach to the 
calculation of the mass function. The reader is referred to Ref. [19] for a more pedagogical introduction. In 
the path integral approach, one treats the variance cr^ = ((5|j) as a "time" parameter, t = i7|j, and considers 
the random walk followed by the smoothed density field 5r as this "time" is increased in discrete steps starting 
from small values (equivalently, as R is decreased from very large values). Here 5r{x) is a stochastic quantity 
in real space due to the stochasticity inherent in the initial conditions. We use the notation Sn to distinguish 



the stochastic variable from the values it takes, which will be noted by Si below. We probe this stochasticity 
by changing the smoothing scale at a fixed location x ^ 0, thus making the variable perform a random walk, 
which obeys a Langevin equation 

!-• (-' 

with a stochastic noise rj whose statistical properties depend on the choice of filter used. In particular, for a 
top hat filter in /c-space, the noise is white, i.e. its 2-point function is a Dirac delta [15], 

{fl{ti)v{t2))^Su{h-t2) (14) 

The random walk can be described as a trajectory {do, 6i, . . . , (5„} which starts with S taking the value Sq = 
at i = (or i? — > oo which is the homogeneous limit), then taking values Si at times i^, finally arriving at 5„ 
at time t„, with a discrete timestcp At = tfc+i ^ tk = tn/n. The probability V{t) that the trajectory crosses 
the barrier at 5c at a time larger than some t (i.e. at scales smaller than the corresponding R or M), is the 
same as the probability that the trajectory did not cross the barrier at any time smaller than t, so that 

r{t)= / d5i...A5nW{{5j}-t), (15) 

J — oo 

where the probabihty density over the space of trajectories, W{{5j}]t) is defined as 

W{{5,}; t) = ( fc(<5(ii) - <5i) . . . 5^{5{U,) - (5„) ) , (16) 

where (5d is the Dirac delta distribution. The first crossing rate is given by the negative time derivative oiV, 
T = —dtV, and the mass function is then / = 2tF{t). In Eqn. (16) one can write the Dirac deltas using the 
integral representation ^(a;) = /_ d\e~'^^^ JI-k, to obtain 

The object ( e^* -^j ^^ ) is the exponential of the generating functional of the connected Green's functions, 
and can be shown to reduce to [37] 
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P=2 ^ ji,..,jp = l 



(18) 



where ( Sj-^ ■ ■ ■ Sj )c is the connected p-point function of 5, with the short-hand notation 5j = S{tj). 

3.1 Halo mass function: Gaussian case, sharp-/;; filter 

In the Gaussian case, all connected n-point correlators vanish except for n — 2, and in the Markovian (sharp-Zc 
filter) case which we are considering, the 2-point function becomes {6j6k) = min(tj,ifc), where imii{tj,tk) is 
the minimum of tj and tk ■ The resulting n-dimensional Gaussian integral can be handled in a straightforward 
way to obtain 

n-l 

M/s- = H *At(4+i - 4) ; ^Atix) = (2^Af)-i/2e-V(2At) ^ (19) 

fc=0 

where we follow MR's notation and use the superscript "gm" to denote "Gaussian Markovian" . As MR have 
shown [19], the resulting expression for 7'gauss(i) in the continuum limit At — > is simply 

Pgauss = 1° d6i... dSr^W^^' = erf ("-^ V (20) 

where we use the notation ly = 5c/(J. (This in principle also includes the redshift dependence of the collapse 
threshold Sc, see below.) This expression for the continuum limit probability T^gauss is of course a well-known 



result going back to Chandrasekhar [38]. This leads to the standard excursion set result for the Gaussian 
mass function /ps = -2i9t|5^7'gauss, 

/psM-Zf^e-"^'/', (21) 

where we use the subscript PS (for Press-Schechter), to conform with the conventional notation for this object. 

3.2 Halo mass function: non-Gaussian case, sharp-A; filter 

In the non-Gaussian case (but still retaining the sharp-fc filter), the probability density W{{Sj};t) also gets 
contributions from connected n-point correlators with n > 3, since these in general do not vanish. These can 
be handled by using the relation A/te* ^j ^ j — —idke^ ^^ ^ % with dj = d/d5j. A straightforward calculation 
then shows the mass function to be 
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where it is understood that one takes the continuum limit At — > before computing the overall derivative 
with respect to t. We will find it useful to change variables from {Sc, t) to {u, t), in which case the partial 
derivative becomes 

-2t{d/dt)\s^ = iy{d/diy)\t - 2t{d/dt)l = i/o>^ - 2tdt . (23) 

It is also useful at this stage to take a small detour and introduce some notation which we will use throughout 
the rest of the paper. We define the scale dependent "equal time" functions 

( S" ) 
e„_2^i^ ; n>3, (24) 

which as we will see, remain approximately constant over the scales of interest. We assume the ordering 
e„_2 ^ ©(e"^^) with e ^ 1, which can be motivated from their origin in inflationary physics, where one finds 
£i ~ /nl^^^^, £2 ~ 9nlA, etc^. Typically we expect e < 10~^ for /nl ^ 100. Fig. 1 shows the behaviour of 
El and £2 in the local and equilateral models, as a function of t = (t|j. We see e.g. that £2 in the local model 
is comparable to e^. In the literature one usually encounters the reduced cumulants 5„, which are related to 
the £n-2 by £1 — crSs, £2 — cr'^S^ and so on. The motivation for using the 5„ comes from the study of NG 
induced by nonlinear gravitational effects. However, as we see from Fig. 1, when studying primordial NG it 
is more meaningful to consider the £„ which are approximately scale-independent and perturbatively ordered. 

We will soon see that a natural expansion parameter that arises in the calculation has the form ^ eiy, and 
we therefore require that the mass scales under scrutiny are not large enough to spoil the relation ev <^\. It 
turns out that observationally interesting mass scales can nevertheless be large enough to satisfy ev^ ~ C'(l). 
Fig. 2 shows the behaviour of £ii^^ and £iv at different redshifts, as a function of mass, in our reference ACDM 
model for local type NG, with f^^^ = 100. The behaviour for the equilateral NG is qualitatively similar. The 
redshift dependence of these quantities comes from the definition of v, 

,(M,z)^V^4^:^^^, (25) 

^'>^ cr(M) D{z) a{M) ' ^ ' 

where we denote the usual spherical collapse threshold as 5co ~ (3/5)(37r/2)^/'^ ~ 1.686, reserving 6c for the 
full, redshift dependent quantity, and a is a parameter accounting for deviations from the simplest collapse 
model. In the standard spherical collapse picture we have a = 1. A value of a different from unity (specifically 



^Notationally we distinguish the order parameter e from the specific NG functions ei and £2. 
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Figure 1: Scale dependence of the £„. Panel (a) : Behaviour of ei vs. a^ in the local and equilateral models, for 
/nl = 100 in each case. Panel (b) : Behaviour of 2nd order [^ e^) terms. We show £2 for the local model with 
/nl = 100 and ^nl ~ 10^. The terms proportional to /^l and (;nl are shown separately. Also shown is Si for the 
same model. 

y/a ~ 0.89) can be motivated by allowing the collapse threshold to vary stochastically [20], as we will discuss 
in section 5. We will soon see that the object ei/^ appears naturally in the calculation, and to be definite we 
will assume ev^ '^ 0(1) for now. In section 4 we will discuss the effects of relaxing this condition and probing 
smaller length scales. 

We now turn to the "unequal time" correlators appearing in Eqn. (22). Since we are concerned with 
large scales, we are in the small t limit, and following MR we expand the n-point correlators around the "final 
time" t. We can define the Taylor coefficients 
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and then expand 
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For the 4-point function we will have an analogous expression involving coefficients Q\ 

Since calculations involving a general set of coefficients C/3, CJ4, etc. are algebraically rather involved, we 
find it useful to first consider an example in which these coefficients take simple forms. In this toy model we 
assume that the e„ are exactly constant, and moreover that the n-point correlators take the form" 



{5f5k5i)c^ei{tjtktiY/'^ ; {5f5k5i5m)c^ e-iitjiktit^f^ . 



(28) 



For clarity, we will display details of the calculation only for this model. We have relegated most of the 
technical details of our calculation to Appendix A. In Appendix A.l wc show that the mass function for this 
model can be brought to the form 
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^Throughout the paper we will consider at most 4-point correlators. This truncation is justified given our assump- 
tions, as we will see later. 
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Figure 2: Panel (a) : Behaviour of e^v"'' 1% vs. mass in the local non-Gaussian model, for /nl = 100. The three 
curves correspond to different redshifts. The horizontal line corresponds to e\v^ j^ = 1. Panel (b) : Behaviour of e\V 
with the same setup as in panel (a). 



where we ignore terms like £\vO(v^'^). The remaining exponentiated derivative can be handled using a saddle 
point approximation. We show how to do this in Appendix A. 2. The expression for the mass function j{y) 
works out to 
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which superficially at least, is comprised of two series expansions, one in the exponential and one as a 
polynomial, both based on the small parameter €v (see however the next section) . 

This derivation assumed that ei and £2 are constant, and that the unequal time correlations are given 
by Eqn. (28). It is straightforward to relax these assumptions and perform the calculation with the exact 
structure of the correlations. Appendix A. 3 shows how to do this, and the result is 
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where the functions c„(i) are defined in Eqn. (A. 17) and characterize the behaviour of the unequal time 
correlations. In our toy model above, the c„ reduce to unity. Indeed, as a check we see that the expression in 
(31) reduces to Eqn. (30) if we take £1, £2 to be constant and set the c„ to unity. 

One issue which we have ignored so far, is that the definition of v involves the variance t = cr^ of the 
non-Gaussian field. Computationally it is more convenient to work with the variance cr^ of the Gaussian 
field in terms of which cosmological NG are typically defined. We should then ask whether this difference 
will require changes in our expressions for /. We start by noting that this difference in variances is of order 
~ e^. For example, in the local model one has cr^{R) — Adi{R) + A{Af^jJd2{R) where A ^ 10^^ is an overall 
normalization constant, di and d2 are scale dependent functions of comparable magnitude on all relevant 
scales, and e is estimated as e '^ JnlA^^^. We therefore have v = 8^1 a = (Jc/cg)(l + ©(e^)). However, with 
our assumption that tv^ ^ ^(1)5 we see that this correction is actually of order ^ {e^v'^)v^'^ ~ e^v^ , which 
we have been consistently ignoring. We will see that even when we relax the assumption ev^ ^ C'(l) and 



probe smaller scales where ev'^ <C 1, this correction can still be consistently ignored. Hence we can safely set 
v ~ Sc/(Jg in all of our expressions. 



4 Consistency of the truncation 

4.1 Comparative sizes of terms in the mass function 

Now that all the derivative operators which we consider important have been accounted for, we can check 
whether our final result is consistently truncated, i.e. whether we have retained all terms at any given order 
in the expansion. Symbolically, our current result for the mass function can be written as 
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-i>.^(l+£!.+£^!.^+0(e^e3,.3)) 
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(32) 



with the understanding that coefficients are computed (but not displayed) for all terms except those indicated 
by the 0() symbols. Also, e^ refers to both e\ and £2- 

Since the expansions involve two parameters, ev and v^"^ ^ they make sense only if we additionally prescribe 
a relation between these parameters. So far we assumed that e is fixed and v is such that tv^ ~ 1, which was 
based on the observation that the term a/' naturally appears in the exponent and is not restricted in principle 
to small values. In Appendix B.l we discuss this condition in more detail, and also analyse the consequences 
of relaxing this condition and probing smaller mass scales. We find that for observationally accessible mass 



scales larger than the scale where tv 



,-3 



, the single expression 
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is parametrically consistent as it stands - the terms ignored are smaller than the smallest terms retained - and 
in fact it remains a very good approximation even when tv^ ~ 1, since the only "inconsistent" term then is 



tv ^, whose effect reduces as i/ increases. On scales where eiy'^ 



and lower, the theoretical error becomes 



comparable to or larger than the quadratic term in the exponential. Plugging back all the coefficients, we 
have the following result for the mass function (excluding filter effects, see section 5), 



.f{i',t) =./ps(i^) exp 
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(34) 



4.2 Comparing with previous work 

In this subsection we compare our results with previous work on the non-Gaussian mass function. As men- 
tioned in the introduction, this quantity has been computed by several authors in different ways [16, 17, 21]. 
If one considers the range of validity of the perturbative expansion, the strongest result so far has been due 

to MVJ [16], who explicitly retain the exponential dependence on £1. 
3 

/mvj = fpsii^)- — rT7 I 1 - -eiJ^ 



as 



1 



Their expression for / can be written 
2rfln£i~ 



3 d\nt 



(35) 



(1 - £i;./3)^/' 

The major shortcoming of their result is that it is based on a Press-Schechtcr like prescription, and must 
therefore be normalized by an appropriate Gaussian mass function, typically taken to be the one due to 



The analysis presented by MVJ in fact allows one to retain terms like 



in the exponential as well, and we 



have seen that when ev ~ 1, these terms are as important as the polynomial eu term retained by MVJ. However, 
since the MVJ expression misses unequal time effects of order ~ ev anyway, it is reasonable to compare our results 
with the expression (35), which is also the one used by most other authors (see e.g. Refs. [39, 40]). 
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Sheth & Tormcn [18]. Additionally, it always misses the contributions due to the unequal time correlators, 
which contribute to the terms ^ ei^, eu^^, etc. in Eqn. (34). When one considers formal correctness on 
the other hand, MR have presented a result based on explicit path integrals, which accounts for the unequal 
time contributions, and which also does not need any ad hoc normalizations (in this context see also Lam & 
Sheth [25]). Indeed, our calculations in section 3 were based on techniques discussed by MR in Refs. [19, 21]. 
As we discuss below however, the fact that MR do not explicitly retain the exponential dependence of eiz^^, 
means that their result is subject to significant constraints on the range of its validity. Their expression for 
/, ignoring filter effects, is^ 

1 3 fi 3 ,^ , 3 / 1 



/MR = /ps(^) ( 1 + TTSi^^' ^ - — (4 - ci) - :^ ( ci - -C2 - 2 ) ^ ) . (36) 



This expression is precisely what one obtains by linearizing our expression (34) in £i, which serves as a check 
on our calculations. LMSV [17] present a result based on an Edgeworth expansion of the type encountered 
when studying NG generated by nonlinear gravitational effects [41]. The result most often quoted in the 
literature is their expression linear in £i (and hence in eiv^), which is 

1 3 f^ 1 /^r, ^filnei^ 2 din El 



/.Msv,nn^/ps(.)(l+ger.-^ 1-^3-2-^ --,-^M- (37) 



In Appendix B.3 of Ref. [17], LMSV also give an expression involving sf and e2, which can be written as 

/LMSV,quad ~ ./ps(i^) 



1 / , , 2dlnei , , 



6 \ ' ly dint 



72^^ ^«(^) + ^Tl^^^^'^) + 24^^ ^^(^) + -.In^''^^^' 



(38) 



where the Hn{i') are the Hermite polynomials of order n. This expression was used by LMSV only as a check 
on the validity of their linear expression. By comparing with our expression which is non-perturbative in 
eij/"^, we will see below that these quadratic terms in fact significantly improve LMSV's prediction. 

Sticking to the linearized results, we see that the expressions of both MR and LMSV have the symbolic 

form 

/ - e~'''/2 [i + £j^3 _^ gj, ^ 1 ^ j ^ (39) 

L ly J 

where the ellipsis denotes all terms of the type ez^~^, ev'^, etc., as well as all terms containing e^. As we have 
seen, deciding where to truncate the expression for / is not trivial, and using our more detailed expression we 
can ask whether the expression (39) is consistent at all the relevant length scales. Immediately, we see that 
this expression cannot be correct once ei'^ becomes close to unity. In Appendix B.2 we discuss the situation 
for MR and LMSV at smaller mass scales. 

5 Effects of the diffusing barrier and the filter 

In Ref. [20] , VIR showed that the agreement between a Gaussian mass function calculated using the statistics 
of random walks, and mass functions observed in numerical simulations with Gaussian initial conditions, can 
be dramatically improved by allowing the barrier itself to perform a random walk. This approach is motivated 
by the fact that the ignorance of the details of the collapse introduces a scatter in the value of the collapse 
threshold for different virialized objects. The width of this scatter was found by Robertson et al. [42] to be 
a growing function of cr(M), which is consistent with the physical expectation that deviations from spherical 
collapse become relevant at small scales. The barrier can thus be treated (at least on a first approximation) 
as a stochastic variable whose probability density function obeys a Fokker-Planck equation with a diffusion 



■'We have corrected a typo in MR's result [21] : the object they define as Va should appear with an overall positive 
coefficient in their Eqns. (85), (87) and (92). 
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coefficient Db, which can be estimated numerically in a given A^-body simulation. In particular, MR found 
Db — 0.25 using the simulations of Ref. [42]. 

Conceptually, the variation of the value of the barrier is due to two types of effects, one intrinsically 
physical and one more inherent to the way in which one interprets the results of simulations. From a physical 
point of view, the dispersion accounts for deviations from the simple model of spherical collapse, for instance 
the effects of ellipsoidal collapse, baryonic physics, etc. On the other hand, the details of the distribution of 
the barrier (and therefore the precise value of Db) will depend on the halo finder algorithm used to identify 
halos in a particular simulation, since different halo finders identify collapsed objects with different properties. 
MR concluded that the final effect of this barrier diffusion on large scales can be accounted for in a simple 
way, by changing i5co — >■ V^Sco where a = (1 + Db)~^- In practice this change is identical to the one proposed 
by Sheth et al. [24]^. As MR argue in section 3.4 of Rcf. [21], this barrier diffusion effect can also be accounted 
for in the non-Gaussian case, again by the simple replacement of Sco — > y/aSco- It is easy to see that their 
arguments go through for all our calculations as well, and we have implemented this change in our definition 
of 1/ in Eqn. (25), setting ^/a = 0.89. 

In Ref. [19], MR also accounted for the non-Markovian effects of the real space top-hat filter, as opposed 
to the sharp-fc filter for which the results of section (3) apply. This is done by writing the 2-point func- 
tion {S{Ri) S{R2)) calculated using the real space top-hat filter, as the Markovian value plus a correction, 
{5{Rj) 5{Rk)) = min(tj, tfe) + A^fe, and noting that the correction Ajk remains small over the interesting 
range of length scales. In fact, MR show that a very good analytical approximation for the symmetric object 
Ajk, is 

A,, ^ Kmin(t„t,) (l Z!!^^') , (40) 

where in our case we find k{R) ~ 0.464 + 0.002i?, with R measured in h~^Mpc. The mass function is then 
obtained by perturbatively expanding in Aij , with the leading effect being due to the integral 



/Sc 1 " 

d5,...d6n- J2 ^jkd.dkW^^, 



J,fe=l 



which on evaluation leads to 

7g,sharp— xt,^i ^J 



2X1/2 



(l-K)e-V2+-r('0,^)+O(«^) 



(41) 



where the subscript stands for Gaussian noise with the top-hat filter in real space, and k introduces a weak 
explicit t{= a^) dependence. In Ref. [21] MR proposed an extension of this result to the non-Gaussian case, 
by assuming that all the non-Gaussian terms that they computed with the sharp-fc filter, would simply get 
rescaled by the factor (1 — k) at the lowest order, but otherwise retain their coefficients. Symbolically, their 
result (Eqn. 88 of Ref. [21]) is 



/ng,sharp— xl,^7 ^j ''^ ^ 



(1 - «)e-^V2 (1 + ,^3 + ^^ + ^^-1) + |r U y ") 



(42) 



with the specific coefficients of the ei'^, ev and tv ^ terms being identical to those in Eqn. (36). However, 
the coefficient of e.g. the nev term arises from the action of an operator ^ ^ ■ j, Ajkdjdk combining with the 
first unequal time operator ^ eit^" X]i(^ ^ ^j)^i Sfc i dkdi, and there is no simple way of predicting its exact 



^A potential issue in this argument lies in the assumption of a linear Langevin equation for the stochastic barrier 
B, resulting in a simple Fokker-Planck equation with a constant Db like the one in MR, while the distribution of B 
was found to be approximately log-normal (and therefore non-Gaussian) in Ref. [42]. One can see that a Langevin 
equation of the type dtB = _B^ (which would produce a log-normal distribution) can be approximated by dtB = (-B)^, 
whenever the fluctuations around (B) are small, and gives a constant diffusion coefficient as long as (B) is constant. 
Although both approximations are reasonable on the scales of interest, non-Gaussian and scale dependent corrections 
to the barrier diffusion should be studied, since in principle they could be of the same order as the other corrections 
retained here. This investigation is left for future work. 
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value beforehand. Since MR explieitly negleet sueh "mixed" terms, their formula is not strictly inconsistent, 
as long as one keeps in mind that the theoretical error in their expression is of the same order as the terms 
~ Ktv that they include. However, if one wants to consistently retain such terms, a detailed calculation is 
needed^. Our calculations (not displayed) indicate that the coefficient of the nev term depends on certain 
details of the continuum limit of the path integral near the barrier, which require a more careful study. We 
are currently investigating methods of computing these effects. At present however, we conclude that the 
mixed terms involving both filter effects and NG, must be truncated at order ^ nev. 

Finally, the filter-corrected mass function is also subject to effects of barrier diffusion. Here we make the 
same assumptions as MR do in Ref. [20], namely that the barrier location satisfies a Langevin equation with 
white noise and diffusion constant Db, which can be accounted for by replacing k ^- k = k/(1 + Db) — an. 
However, it is difficult to theoretically predict the unequal time behaviour of the barrier correlations, and 
these simple assumptions must also be tested, perhaps by suitably comparing with the detailed results of 
Robertson et al. [42] . We leave this for future work. Our final expression for the mass function, corrected for 
effects of the diffusing barrier and the top-hat real space filter, is 

1 3 1/2 e2\ 4 

6 8 V 1 3/ 



/(z.,t)-/psM l-S + 0(^^) exp 



(I — 2dliiK,/dlnt) 9/ Ts 1 / ci 

1 + -5^ '- ^-kiy-^ (l - 2iy-^) - -£iv — — + 4 - 2ci 

1 — K 4 \1 — K 

1 ,_1 / 1 



~-e^v- (^ci--C2-2 

+ 0{k^i^-^, kev, kv-*^) + C'(e2^/^ e3^/^ ev-^) \ , (43) 

where we have chosen to account for the scale independent 0{k'^) error arising from filter effects, as an overall 
normalization uncertainty, and have explicitly displayed the orders of the various terms we ignore. Here 
fpsi^) is given by Eqn. (21) with v{M,z) defined in Eqn. (25). 

To summarize, Eqn. (43) gives an analytical expression for the non-Gaussian mass function. This expres- 
sion is based on approximations that are valid over a larger range of length scales than the ones presented by 
MR and LMSV, and incorporates effects which are ignored in the expression presented by MVJ and LMSV. 
Like all these other mass functions, it suffers from the errors introduced by filter effects. However, the largest 
of these can be accounted for as an overall normalization constant, which can be fixed using, for instance, 
results of a Gaussian simulation. Among the NG functions £i, £2 and the c„ which appear in the mass func- 
tion, the most important ones over the mass range of interest are £1, £2 and ci which are nearly constant. In 
principle though, all these functions must be computed numerically for every mass scale of interest, and indeed 
all the plots in this paper use the results of such numerical calculations. However, since this is somewhat 
tedious to do in practice, in Table 1 we provide analytical approximations for £1, £2, ci, C2 and C3, for the 
local and equilateral case as a function of cr^. As mentioned earlier, all these quantities are independent of 
rcdshift, although they depend on the choice of cosmological parameters in a complicated way in general due 
to the presence of the transfer function in their definitions. However, the dependence on erg is simple, and 
one can check that we have £1 oc erg, £2 oc cr| and that the c„ arc independent of erg. Recall that the c„ 
are also independent of /nl and ^nl- For completeness, in Table 1 we also give approximations for the filter 
parameters k and dlnk/dlnt which appear in the mass function, as functions of a^ . 

6 Results and Discussion 

In this section we conclude with our final results for the non-Gaussian halo mass function, comparing our 
approach with previous work. In principle, we should compare the full expressions for the mass functions 
of various authors with ours. However, recall that for MVJ and LMSV one has to multiply an analytically 



^Notice that this issue is completely decoupled from the subtleties in truncation discussed in section 4 ~ this problem 
remains even at scales where MR's expression is formally consistent. 
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Parameter 


Fitting 


form b + cr 


Local NG 


b 


c 


n 


£1 


0.0096 


0.015 


0.18 


Cl 


0.98 


0.073 


0.094 


C2 


3.15 


0.79 


0.69 


C3 


2.15 


0.45 


0.65 


£2(/Sl) 


-0.0049 


0.0059 


0.011 


e2(5NL) 


7.9 • 10-4 


0.0022 


0.25 



Parameter 


Fitting form & + ci" 


Equilateral NG 


6 


c 


n 


£l 


0.01 


-4-10-4 


1.25 


Cl 


1.03 


-0.052 


0.30 


C2 


2.32 


0.93 


0.49 


C3 


1.72 


0.36 


0.54 


Filter 




k 


0.36 


0.015 


-0.47 


d\nk/d\nt 


0.046 


-0.064 


-0.17 



Table 1: Analytical approximations for the various NG parameters, in the local and equilateral cases, as a function 
of t = cr^, in the range 2 ■ lO" < M/ih'^M^oi) < 5 ■ 10^^, for /nl = 100 and ^nl = lO". We have ei oc /nl in both 
cases, and for £2 in the local case we give separate approximations for the terms proportional to /^l and ^nl. We 
do not consider £2 in the equilateral case, since the trispectrum in this case is highly model dependent. We also give 
approximations for the filter parameters R and d\nk/d\nt as functions of i, in the same mass range. The errors on all 
the approximations are less than 1%, except for e2(/NL) where the error is ~ 6%. This was due to numerical difficulties 
in calculating this object. These approximations of course depend on our choice of cosmological parameters. 



predicted ratio i?ng — f{i>, M, /nl)//('^, M, Inl = 0) with a suitable Gaussian mass function based on fits 
to simulations. It is not clear how to compute theoretical error bars on the latter. On the other hand, the 
object Rng itself is an unambiguous theoretical prediction of every approach, that is MVJ, LMSV, MR and 
our work, and we can compute theoretical errors on it. In this work, we will restrict ourselves to comparing 
the different expressions for R^g- In future work, we hope to compare both _Rng and the full mass function 
with the results of iV-body simulations. 

In Fig. 3 and Fig. 4 we plot the ratio Rng, respectively without and with the filter effects, at redshift 
z = 1. In this way we can explicitly disentangle the errors due to an approximate treatment of non-Gaussian 
effects, from those due to the filter effects. We compare our expression (43) with the expressions of MR (36), 
LMSV (37) and (38), and of MVJ (35). Notice that, when considering the filter effects, the Gaussian function 
that enters in the ratio R^g is defined to be the function with /nl = (i.e. without NG but with filter 
effects when present). We use the local model, setting /nl = 100 and ^nl = 0, and use the reference ACDM 
cosmology described in section 2. We do not explicitly show the final results for the equilateral model, but 
they are qualitatively similar. As is commonly done in the literature, we modify the LMSV and MVJ curves 
by applying the Sheth et al. correction of 5cq -^ y/adco- An identical correction is already present in the 
expressions (43) and (36) due to the barrier diffusion. We set y/a — 0.89, which is the value inferred by MR in 
Ref. [20] using the simulations of Robertson et al. [42] . We wish to emphasize a feature that our calculation 
shares with that of MR, which is that the constant a is the only parameter whose value depends on the output 
of A'^-body simulations. The rest of the calculation for the mass function is completely analytical and from 
first principles. 

To make the comparison meaningful, we introduce theoretical error bars on the curves. These error bars 
have no intrinsic statistical meaning - they simply keep track of the absolute magnitude of the terms that 
are ignored in any given prescription for computing the mass function. As we have discussed at length in 
section 4, these theoretical errors are scale dependent. The estimated error magnitude for each point is the 
maximum among the terms ignored in the expression. More explicitly, the errors for the linearized LMSV 
expression (37) are estimated as the maximum of (eu^)'^ which comes from the expansion of the exponential, 
eiy which is the order of the largest unequal time terms missing, and ku^^ which comes from the filter effects. 
The errors for the LMSV expression (38) are similarly estimated as the maximum of (ei^'^)^, ev and ku^^. 
The largest error for the MVJ expression (35) is the maximum of ev (unequal time terms) and ku~^ (filter 
effects). Finally, the error for the MR expression (36) is the maximum of (eiy^)^ from the expansion of the 



exponential, eiy ^ from the largest unequal time terms ignored, and k^v 
include the filter effects and the associated errors only in Fig. 4. 



and kei' from the filter effects. We 
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Figure 3: Theoretical corrrparison of the different mass functions at 2 = 1, without the filter effects, i.e. setting 
R — 0. We plot the ratio _Rng of the non-Gaussian and Gaussian mass functions, in the local model with /nl ~ 100 
and gfNL ~ 0. See main text for a discussion of the error bars. The arrow indicates the mass scale where eii/^/6 = 1, 
i.e. where the expansions of LMSV (both linearized and quadratic) and MR break down. 

From these figures, we can draw some interesting conclusions. First of all, we see that it is important 
to retain terms which are quadratic in the NG, either with a saddle point method like in MVJ and in our 
formula, or by expanding the exponential up to second order, like in LMSV. Actually, we argue that it is 
correct to keep the exponential, otherwise the expansion breaks down when ev^ is of order unity. We notice 
in passing that the term proportional to £2 which comes from the trispectrum, partially cancels with the e^ 
term. Secondly, comparing our expression with MVJ's, we can observe that keeping the unequal time terms 
allows us to sensibly reduce the theoretical errors due to the approximate treatment of NG. In fact, if these 
terms are missing, they provide the largest theoretical error on large scales. Instead, the largest theoretical 
error on small scales comes from the approximations involved in dealing with a real space top-hat filter, as is 
apparent from Fig. 4. 

To conclude, in this work we have calculated the non-Gaussian halo mass function in the excursion set 
framework, improving over previous calculations. We started from a path integral formulation of the random 
walk of the smoothed density field, following Ref. [19]. This allows us to take into account effects due to 
multi-scale correlations of the smoothed density field ("unequal time" correlations), and due to the real space 
top-hat filter, which generates non-Markovianities in the random walk. We recognize two small parameters 
in which we perturb: e, defined below Eqn. (24), which measures the magnitude of the primordial NG; and 
u~^ = au/Sc, which is small on very large scales. In order to do a consistent expansion and to estimate the 
theoretical errors, one must study the (scale dependent) relation between these two parameters, which we 
have discussed in Sec. 4. We then used saddle point techniques which allowed us to non-perturbatively retain 
the dependence on ev^, which naturally appears in the calculation and whose magnitude becomes of order 
unity at high masses and high redshift. Finally, we included effects due to the choice of filter function and 
due to deviations from spherical collapse, as explained in Sec. 5. Our final result is presented in Eqn. (43), 
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Figure 4: Same as Fig. 3, but including filter effects. Tfrese affect only the error bars for MVJ and LMSV, and they 
affect both the curve and the error bars for MR and our result. For MR and our result, the Gaussian mass function 
used to construct the ratio -Rng, is taken as the non-Gaussian result at /nl = 0, and hence includes filter effects. 



which we reproduce here: 

/(^, = /psM f 1 - « + 0(^') 1 exp 
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In Table 1 we provide analytical approximations for the various parameters that appear in this expression, 
which in general must be computed numerically. We also considered other expressions for the mass function 
found in the literature, which use different expansion methods but do not estimate the theoretical errors. 
We estimated the theoretical errors for each formula, and we show comparative plots in Fig. 3 and Fig. 
4. In our work we have improved over the calculations of MVJ [16] and LMSV [17] (who ignore unequal 
time correlations) and of MR [21] (who do not retain the exponential dependence on €v^). We have also 
demonstrated that the (linearized) result of LMSV can be significantly improved by retaining the quadratic 
terms of their calculation which are usually ignored in the literature. We find that at large scales and high 
redshifts, the biggest theoretical errors are introduced by ignoring the exponential dependence on ev'^ , followed 
by the neglect of unequal time correlations. The errors on our expression (44) are therefore significantly smaller 
than those of the others. The strength of our approach lies in the combination of path integral methods as 
laid out by MR [21], and the saddle point approximation as used by MVJ [16]. 

Our work can be continued in several directions. First, a thorough calculation of the effects due to the 
choice of the filter should be performed, since they lead to significant uncertainties in our final expression. 
This would include a study of the details of the continuum limit of the path integral near the barrier, and also 
a study of the statistics of the barrier diffusion process in the presence of filter effects. Second, a comparison 
with A^-body simulations should be performed, in order to quantitatively assess the possibility of constraining 
NG using our work. Third, it would be interesting to study how to account for the effects of ellipsoidal 
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collapse, in a framewrok such as the one employed in this paper. Finally, an application to the void statistics 
along the same lines should be feasible. The problem here is made more interesting by the presence of two 
barriers, as discussed by Sheth & van de Weygaert [43], and since voids probe larger length/mass scales than 
halos, they constitute a promising future probe of primordial NG [44]. 
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Appendix 

A Mass function calculation 

A.l Equal time vs. unequal time terms 

In this appendix we show how the exponentiated derivative operator in the path integral can be handled by 
separating the contributions of the equal time and unequal time correlations. While this calculation assumes 
the toy model introduced in Eqn. (28), it easily generalizes to the more realistic case as we discuss later. 
Using the first few terms of the unequal time expansions, in our toy model one can write 






j,k,l=l ^ j.,kJ=l 3 = 1 ' kd = l 



SjSkSiSm )c djdkdidm = £2^^ f X! djdkdidm -2^(1- ^)dj ^ dkdid„i + ••■]■ (A.2) 



j.k.l,m=l ^ j,k,l.m=l i— 1 fc,/,r?i— 1 

These derivative operators are exponentiated in the path integral, and act on W^"^. One simplification that 
occurs in our toy model, is that the path integral in Eqn. (22) becomes a function only of v (although this 
is not obvious at this stage), and hence eventually only the i^S^ part of the overall derivative contributes. 
However, the structure of the exponentiated derivatives is still rather formidable. Moreover, the truncation 
of the series at this stage is based more on the intuition that higher order terms should somehow be smaller, 
rather than on a strict identification of the small parameters. In fact, we will see in detail in section 4 that 
the issue of truncation involves several subtleties. 

To make progress, it helps to analyze the eflfect on W^"^ of each of the terms in the above series, before 
exponentiation. The leading term in Eqn. (22) involves the multiple integral of W^"^, which is just the 
quantity T'gauss encountered in Eqn. (20). The operator vdi, acts on the error function to give the Gaussian 
rate of Eqn. (21). Next, notice that the action of the operator X]i=i ^i "-"^ any function g{Si, . . .dn) under 
the multiple integral, is simply 



<5<: "a rS. 



d 
dSi... dSn ^ djg = — / dSi . . . dSng , (A.3) 



j = l ---^-oo 



Using this, and the fact that t^''^{d/dSc)\t — du\t, we see that the leading term in Eqn. (A.l) (i.e. the term 
with no powers of (1 — tj /t)), leads to a term involving 



£1 



vd.{d,fcii (v/V2) ^ fpssiiy^l + 0{iy-^)) , 



The problem with this term is that the quantity eiiy^ can be of order unity, and hence cannot be treated 
perturbatively. To be consistent, we should keep all terms involving powers of Sii^"^. Luckily, this can be 
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done in a straightforward way due to the result in Eqn. (A. 3). We see that the entire exponential operator 
exp[—{eit^^^ /3\)^^ f, i^^djdkdi] in Eqn. (22) can be pulled across the multiple integral and converted to 
exp [—{ei/3l)df,] acting on the remaining integral. Similarly, the operator exp [(e2t^/4!) X]?fc / m.=i djdkdidm] 
can be pulled out and converted to exp [{e2/4:\)d^], and the same applies for all such equal time operators. We 
will see later that the action of these operators can be easily accounted for, using a saddle-point approximation. 
To summarize, the function / at this stage is given by 

j=l k,l=l j-k=l 1=1 

1 / " t " 

-i!^2t'(2E(l-f)^^- E dudid^ + . 

^ j — 1 k,l,m—l 



W^sm (A.4) 



Now consider the action of the individual terms in the remaining exponential under the integrals, but without 
exponentiation. From MR [21], we have the following results^, 

E (1 - f ) E / ^'^1 • • • d6nd,dkdiW^^ = - Ty^e-"'^^ , (A.5a) 

j=i ^ kj=i-'~^ ^'^^ * 

n f '"■ r^<: / 9 \ ^/^ % 

T.(^-i)'T. ^'^1 ■ • • d6nd,dkdiW^^ = - -^h{^) , (A.5b) 

j=l kj=l''~°^ ^""^ 

" , /, " /"^c /9\^/^ 4 

E(l-^)(l-y)E/ d5,...d5r.d,dAW^-^i-) ^h{.), 

i.k=l J=l "'-°o ^""^ 



(A.5c) 
where we have defined 

where r(— 1/2, i^^/2) is an incomplete gamma function. Let us focus on the term in Eqn. (A. 5a). If we 
linearize in ei in Eqn. (A. 4), then this term appears with eii^'^S^ acting on it, leading to ~ fps£iv <C /ps- 
This term can therefore be treated perturbatively. Similarly, one can check that the terms given by Eqns. 
(A. 5b) and (A. 5c) also lead to perturbatively small quantities, which are in fact further suppressed compared 
to £iv by powers of v^^ . Specifically, one obtains terms involving eicrfc (j^/v2) which, for large v, reduces 

A few comments are in order at this stage. First, this ordering in powers of v^^ is a generic feature of 
integrals involving an increasing number of powers of (1 — tj /t) being summed. This can be understood in 
a simple way from the asymptotic properties of the incomplete gamma function, as we show in Appendix C. 
We are therefore justified in truncating the Taylor expansion of the unequal time correlators, even though 
superficially (on dimensional grounds) each term in the series appears to be equally important. Secondly, we 
have not yet accounted for the effect of the exponential derivatives. In fact we will see in the next section 
that when ev^ ^ ^(1): it is these terms that impose stricter conditions on the series truncations. For now, 
however, we have no guidance other than the fact that if we account for one term of order ^ e^v"' , then we 
should account for all terms at this order. Given this, note that for eu^ ^ C'(l) we have ly^^ ^ eu, and hence 
the terms arising from Eqns. (A. 5b) and (A. 5c) are of order ^ e^v'^. To consistently retain them, we must 
therefore also retain the term linear in £2 and the one quadratic in ei, when expanding the exponential. These 



^The terms in Eqns. (A. 5a), (A. 5b) and (A. 5c) are, upto prefactors, the integrals of what MR denote as n'"^''^'"-', 
n(3,NNLa) ^^j n'^'NNLb) respectively in Ref. [21]. 
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involve the following quantities: 

r^c /o\i/2 



^(1-f) Y. / d5^...d5ndjdkdidraW^^ = -[-] -^^ e'^ /\ (A.7a) 

3 = 1 * k,l.,m=l-'-°° ^^^ * 

^(i-|)(i-^) ^ / dS,...d5^d,dkdidiAAsW^'^ = -[-) Tj^e--^'/^, (A.7b) 

where we have used the result (A. 3), and in Eqn. (A. 7b) also the identity 

dlh{iy) = -iye~'''/^. (A.8) 

We now see that the result of the path integral depends only on i/. Putting things together and computing 
the overall i^ derivative, we find the result in Eqn. (29). 

A. 2 Saddle point calculation 

To compute the action of the exponentiated derivative operators, we start by writing the expression in square 
brackets in Eqn. (29) in terms of its Fourier transform, using the relations^ 

J-co vStT 
J-ca V 27r 
7_no V 27r 



1'\^/^ - ''ji_\= /"°° _^l„»Ai. -AV2 



crfc ^ = / ^-e'^'^e-^ '' . (A.9) 

2^ V\/2y i-ooV2^A ^ ' 

Together with the identity ^J^i-^ld'^T ^i^'^ ^ ^A(-i\)^ ^i\v ^ j^j. constant A and B, this gives 

/(^) =(iy\ r 4^e.A.e-AV2+(-.A)3si/6+(-^A)^e./24+...p(;^) (^.10) 

XT^J J-oc V27r 



where 'P(A) is the truncated series given by 

The integral in eq. (A. 10) can be performed using the saddle point approximation. We write it as 

f^,)^(lY\ r l^e'^W, (A.12) 



."^y J-oc V27r 

where 

0(A) ^ zAz/ - Ix^ + ^A3 + £^a4 + ln7'(A) + . . . (A.13) 

2 D 24 

The location of the saddle point, A = A*, is the solution of (f>'{X^) = 0, and the saddle point approximation 
then tells us that 

^^=e*(^) = e^(^*)(|c/."(A*)|)-i/2 ^ (A.14) 



2tt 



*We are using a regulator which shifts the pole at A = in the last expression in Eqn. (A.9), to A = —ia where a 
is real, positive and small. 
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(see Appendix D for a discussion of the errors introduced by this approximation). It turns out that in order 
to obtain f{v) correctly up to order ^ e^u^, we only need A* correct up to order ^ eiy. The expression for 0' 
at the relevant order is, 



(/.'(A) = 21/ - A + ^A^ 
and solving for A, perturbatively up to order ev, we find 



A, = w 
This leads to the expression for /{u) in Eqn. (30) 



l-^eiiy + 0{e^iy^) 



(A.15) 



(A.16) 



A. 3 Result with full unequal time terms 

In the more realistic case of slowly- varying e„, we choose to parametrize the coefficients Q3 and Qi (see Eqn. 
(26)) in a convenient way as follows : 

G!,'^'-'' - U{t)c^it)t^/' ; g^s"^'-'' = -isi(t)c,(t)t-V2 , 



i'^'-'^ = ]e^{t)c,{t)t-^/^ ; g^^-'-'-'^ = U{t)c,{t)t, 



(A.17) 



where the coefficients c„(i) are smoothly varying functions and depend on the NG model. They are defined 
in such a way that they all reduce to unity in the toy model defined by Eqn. (28). Fig. 5 shows the behaviour 
of ci, C2 and C3 with a^ , for the local and equilateral models. The £„ and c„ are independent of redshift by 
construction, since the linear growth rate D{z) always drops out in their definitions. Further, the c„ do not 
depend on the values of /nl and .9nl- One can then use the definitions of £1 and the c„ to prove the following 
useful relations 



dlnei 3 dlnci 

d\nt ~ 2 ^^^ ~ -* ' d\nt 



_, 3 1/1 



(A.18) 



The calculation of the mass function for this general case proceeds completely analogously to that for the toy 
model, apart from a few subtleties which we will discuss later. In this case Eqn. (A. 4) is replaced with 



/ = (/.9. - 2ia0e-(^^(*'/^''^'^+(^^(*)/^')^'+-g(i/,i) 



1 dhiei 2 1 din 82 ^ 3 



3 dint 



12 dint 



2ig-(e.(t)/3.)a^ + (e.(t)/4.)8^ + ...3^^(^^^)^ 



(A.19) 



wher the function g{i',t) can be shown to be 



gi^.i) 



1/2 



7r\i/2 / V 



erf 



VV2 



1 



-eicie 



-v^/2 



£1 /3 



-C2 - 2c3 h{iy) 



\e\c\ve-^'l^ + ^e2C4^e-V2 



(A.20) 



The expression in Eqn. (A.19) can be evaluated analogously to Eqn. (29), since the additional derivatives pose 
no conceptual difficulty. The result of the saddle point calculation, correct up to quadratic order assuming 
tv^ ~ C(l), and after using the relations (A.18), is given in Eqn. (31). 
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Figure 5: The derivative coefficients ci (panel (a)), C2 (panel (b)) and ca (panel (c)), as a function of a , for local 
and equilateral NG models. These quantities are independent of redshift and the NG amplitudes /nl and jtnl- The 
axes are logscale. 

B Truncation of the perturbative series 



B.l Analysis of "transition points" 

In this appendix we analyse the consistency of our series truncation at various mass scales. When eiy^ 



1, 



the polynomial in (32) we retain the terms ep c:i u , {ev , e v ) 



, and we discard (ei^ ^,e ,e^iy^) 



m 

'6 



It would seem that our expression is then correct upto order ~ iy~ . However, the terms discarded in the 
exponential have the form exp(0(e'^j^^)) ~ exp{0{i'^*)) ^ 1 + 0{v~^). The error we are making is thus of the 
same order as the smallest terms we are retaining, and it therefore makes sense to also ignore all the terms of 



order 



which we computed in the polynomial. The consistent expression when ev ~ 1 is then given by 



/ ^ e-5-'(i+-+-'-') [l + eu + O {v-^)] 



(B.l) 
"2, etc. It 



Clearly, similar arguments can be applied at smaller scales where, e.g. one might have ev^ ~ v^ 

is then important to ask which mass scales correspond to these "transition points" . In Fig. 6 we plot v{M^ z) 

given by Eqn. (25) in an observationally interesting mass range, for three different redshifts. The horizontal 



-1 



V 



V 



V 



lines mark the transition points where (.v^ becomes equal to (from top to bottom) 1, v 
and v~^ . We fix e = 1/300 which follows from the fact that in the local model with /nl = 100 we have 
£1 ~ 0.02 (see Fig. 1), and the expression for J(y,M) contains the quantity ei/6 in the exponential. From 
the intersections of the horizontal lines with the curves, we see that different transition points are relevant 
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Figure 6: i/ = 5c(2)/o-(M) in the range 5 ■ lO" < (M//i"^Msoi) < 5 • lO^'^ for three diflterent redshifts, with e = 1/300. 
The sohd, long dashed and short dashed curves correspond to redshifts 2 = 1, 0.5 and respectively. The horizontal 
lines mark the transition points where ev'' becomes equal to (from top to bottom) 1, v"^ , v~^ , i^~^, i^~* and u~^. 

at different redshifts, and their locations also obviously depend on the value of e. For example, we find that 
the transition point where ev^ ~ i'"^, remains accessible even when e is an order of magnitude smaller (with 
e ~ 1/3000, this transition occurs at i/ ~ 4.96). The transitions at ev^ ~ 1, i^^^ on the other hand, are not 
accessible for this level of NG. The transition at eiy^ ~ ly^^ is therefore observationally very interesting. 

We will now discuss in some detail the truncation of our expression for /, at various transition points. 
The goal is to try and settle on a single expression which is valid over a wide range of scales (i.e. across 
several transition points). This can then be applied without worrying about truncation inconsistencies. Of 
course, the order of the discarded terms will then depend on the particular transition point being considered, 
leading to a scale dependent theoretical error. 

B.1.1 ei/3~i/-i 

At this transition point, the terms we retain in the exponential are 

while discarding 0{e^iy^) = 0{h'^ ). In the polynomial meanwhile, we retain 



-3 



£2^.2 ^ ^-6 



while discarding 



0{€iy-') = 0{iy-') ■ O(e^) = 0(z/-») ; 0(e^i/^) = 0(z/-^) 



Our expression (32) therefore retains all terms correctly up to order ^ v ^ , and is consistent. With some 
foresight however, it turns out to be more convenient to degrade this expression somewhat by also discarding 



the polynomial quadratic term e v 



2,, 2 



The remaining expression, 



/ ^ ^-k-'{^+--+-'-') fi + e^ + 1 + o (^-6) 



fB.2) 



is also consistent at this transition point, and has a form which is identical to the ones we will see next. 



B.1.2 eu^ ~ u-^ 

As we mentioned earlier, this transition point is observationally quite interesting. The terms we retain in the 
exponential are 



eu^ ~ v-^ 



e^u^ 
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while discarding 0{e ly'^) = 0{v ), and in the polynomial we retain 



v-^ ■ tv-^^v-'' ■ e2j.2~j.-8. 



while discarding 

This time we see that the term ev^^ has become as important as the quadratic term f^v^ in the polynomial, 
and to be consistent we should discard the quadratic term. The expansion should read 



/ ~ e-^-^(i+-+^^-^) [i + ei. + 1 + (^,-8)1 . (B.3) 

L V \ 

B.1.3 ei/3 ~ jy-3 

A similar analysis as above shows that at this stage tv^^ ~ v^'^ > e^j/^^ q^^A a consistent expression again 
requires dropping the quadratic term in the polynomial, leaving 

/ ~ e-^^(i+-+e^-^) [i + et. + 1 + (^,-9)] (B.4) 

B.1.4 ev^ ~ i/~^ and smaller 

Beyond this point, the term ev^^ which we discard in the polynomial, becomes comparable or larger than the 
quadratic term of the exponential as well, and a consistent expression becomes 

/~e-5'^'(i+<^'^) [l + £!/+- + ... 1 (B.5) 

L V J 

The parametric order of the terms now discarded, depends on the exact relation between ev^ and v~^ . 

Finally, note that the error introduced by setting v -^ v^ where v^ is defined using the variance of a 
Gaussian field, was estimated in section 3 as 0{e'^). When ev^ ~ 1, this error is of order 0{e^i'^) and can 
therefore be consistently ignored. It is not hard to see that at all lower transition points, this error continues 
to be comparable to or smaller than the largest terms being discarded, and can hence be consistently ignored. 
This finally leads to the conclusion stated in the main text. 

B.2 Truncation in MR and LMSV results 

Since the mass scale where ev^ ~ 1 is on the border of the observed mass window (for galaxy cluster ob- 
servations), even at high redshifts, let us therefore directly look at the case eiy^ ~ i^^^ which, as we saw, is 
accessible over a wide range of redshifts for e ~ 10"^^ and at high redshifts also for e ~ 10^'^. In this case the 
terms MR and LMSV retain have magnitudes 

and terms like ei^^^ ~ j/^* are discarded. We know from our expression however, that ei/^ appears in 
the exponential, and therefore leads to terms like {ei'^)'^ ~ ly^^ and (eiy^)^ ~ i/~^ when the exponential is 
expanded, which are of the same order as the terms retained in (39). The exponential also contributes a term 
g2j^4 ^ ^-6^ which in fact involves the irzspectrum of NG, again at the order retained by MR and LMSV. 
The error in the expression (39) when ev^ ~ i^~'^, is therefore 0{€v). (A similar analysis shows that the error 
at transition point where ev^ ~ v^^ ^ is 0{v^'^) > 0{eiy).) 

From a purely parametric point of view, the situation for MR and LMSV improves as v is decreased 
further, and the expression (39) as it stands, becomes exactly consistent (in the sense discussed in the previous 
subsection, see below Eqn. (33)) when ei^^ ~ i^~^, because at this stage ei'~^ ~ i'~^ while {eiy^)'^ ~ zy^^^ and 
g2j^4 ^ zy~i2^ ^jj(j hence the exponential only contributes a single linear term ei^^. More importantly, LMSV's 
expression also has errors due to the absence of the unequal time terms discussed earlier, which are of order 
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^ ev and can be dominant over the others. For the intermediate transitions, the analysis shows that when 
gjy3 ^ j^-3^ ^j^g error in (39) is 0{v^^) > 0{€v^^), and when ev^ ~ v~^ , the error is 0{ev~'^). This should be 
compared with our result (34), in which the error (at least on large scales) is always parametrically smaller 
than the smallest terms we retain. 

C Hierarchy of terms in Eqn. (A. 4) 

Here we argue why the hierarchy of terms ordered by powers of v^"^ emerges on expanding the exponentiated 
derivative operators in Eqn. (A. 4). Focusing on terms involving the 3-point correlator, one sees that a generic 
term in the expansion contains some powers of (eii'^'^), multiplying an n-dimensional integral containing 
some summations ~ ^" . =i(l ^ tji /t)^^ (1 ~ ^j2/0^^ • • • ^ji^j2 • • •: ^'^'^ ^l^o some summations over "free" 
derivatives ~ ^^ j. ^-^ dkidk2 ■ ■ •, all of this acting on W^™. More precisely, the structure of the terms is 

X [(1 - t,,„+, /t) ... (1 - i,3™ /t)] ^ d,, . . . %,„ W^«- , (C.l) 

for TO > 1 and non-negative p, q, r such that not all three are zero. The terms we have considered in the text 
are {m,p,q,r) = (1,1,0,0), (1,1,1,0), (1,2,0,0) and (2,1,0,0). We have already discussed how the "free" 
derivatives can be pulled out of the integral and converted to d^. For the "non-free" derivatives, we see that 
what is important is the total number of (1 — tj/t) factors accompanying these derivatives. For example, 
the (1, 1, 1, 0) term in Eqn. (A. 5c) has the same structure as the (1, 2, 0, 0) term in Eqn. (A. 5b) - the effect 
of J2i fc(l ^ ^i/^)(l ^ tk/t)djdk, up to numerical factors, is identical to that of ^ ^,(1 — tj/t)'^djdk- This is 
expected to be true also with higher numbers of non-free derivatives. 

It is then possible to understand the hierarchy of terms by only considering terms containing "^Al — 
tj/t^dj, and no other non-free derivatives. The basic object to study now becomes 



T.(^-t,/t)J 



d5i...d6„d,W^" 



which in the continuum limit can be shown to reduce to the integral 

Notice the similarity with the function h{i') in Eqn. (A. 8), which of course is not accidental given the 
definitions of these objects. It is now easy to check that increasing the powers of (1 — tj/t) in the summation 
amounts to increasing the powers of (1 — y) in g/Qy We are then comparing (with A = v'^/2) g/Q\{A) with 
g{p){A) where 

miA)^f^^,{l-yf'^^^e-^'^. (C.3) 

Starting with p ~ I and manipulating the integrals, it is straightforward to establish the recurrence 



5(p+i) (A) = g(p) (A) - / dA g^p) [A) . (C.4) 

J A 

The argument is now almost complete. We know that for large A ~ v'^ l'2,^ we have T{n^A) ~ er^A^^^{l + 
0{A^^)). Hence g{Q){A) = {^/2)A~^/'^e^^{l + 0{A^^)), and its integral from A to oo gives a leading term 
proportional to r(-3/2, A) = e-'^A-^l'^{l + 0{A~^)). The pattern is now clear: 5(p)(A) - A-3/2-Pe-^(l -F 
0{A^^)), and since A — i^^/2, this explains the hierarchy of terms in powers of t^^^, in Eqn. (A. 4). 
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D The saddle point approximation 

In this appendix we discuss the saddle point approximation of the integrals of the type appearing in section 3.1, 
and estimate the error it induces. We will argue that the errors introduced by the saddle point approximation 
are much smaller than those due to truncating the perturbative series in the small parameters e and ly^^. For 
an introduction to the saddle point approximation see Ref. [45] . Since we only wish to discuss the saddle point 
method in this appendix, we will ignore here the complications introduced by the unequal time correlators, 
i.e. in Eqn. (A. 10) we set V{\) ~ 1. We will also work here to first order in ev. The extension to a more 
general case is straightforward and the result is given by (31) as described in section 3.1. We begin with 
expression (A. 10): 

/(.) = f^V'. r ^e^W , (D.l) 



where g{X) = ivX - \^ /2 + {-i\fei/Q + ©(e^A*). 

We first find the location of a saddle point A* of the function g{\), by perturbatively solving g'{\*) = 
using €v as the small parameter and demanding g"{X*) < 0. The first-order solution is 

A, = ijy(l - eii//2 + ©(eV^)) , (D.2) 

.9(A*) = -y(l-^ei^ + 0(6'^')), (D.3) 

g"iX,) = -l-eiiy + 0{e^iy^) . (D.4) 

The saddle point approximation consists roughly of performing a Taylor expansion of 5(A) to second order 
around A* in the integrand of (D.l) and performing the resulting Gaussian integral. We will carry this out 
explicitly below. The saddle point prescription will give a good approximation to the integral as long as g{X) 
attains a global maximum at A* (along the contour of integration) ; this is indeed our case since the integrand 
in Eqn. (D.l) will be nearly a Gaussian centered at A* in the complex plane. 

Notice that Im A* ^ 0, requiring a deformation of the contour of integration such that it passes through 
A*. The deformation of the path of integration can be performed by taking a closed contour formed by four 
pieces: The real axis Ci, the line Im A = Im A* which we call here — C2, and the closures of this contour at 
possitive and negative infinity. The integral in this closed contour must be zero, and since the integral on the 
closures of the contour at infinity can be assumed to vanish, we have J^ ~ In ' Therefore C2 is the desired 
deformation of the contour which passes through A*^. We can then make a series of approximations in the 
integral (D.l), which we discuss below, 

,S(A.)+3"(A.)(A-A.)V2 
,-1/2 
^ g-i.^(l-ei./3+0(e^.^)) (^ ^ ^^^ ^ 0(e2,.2))-l/2 _ ^^_^^ 

Here the integrations are performed along the deformed contour. 

In order to estimate the errors induced by the approximation done in equation (D.5), one can keep higher 



V2^ 


i-00 V27r 


:: 


= e«(^*) (-.9"(A,)) 



^Technically, one should also require that Im g{X) be nearly constant along the deformed contour for the saddle 
point approximation to work. In our case one can show that Im g will be suppressed by e. This and all errrors induced 
by the saddle point are accounted for in equation (D.6). 
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Figure 7: Panel (a) : Fractional difference between the saddle point approximation on the r.h.s. of Eqn. (D.7) and 
the numerical integration of the l.h.s. of the same equation. Panel (b) : Total error induced on the result of the toy 
model (D.7) by both the saddle point approximation and the perturbative expansion to leading order in ev. We plot 
the fractional difference between the numerical integration of the l.h.s. of Eqn. (D.7) and the approximation (D.5). 
Both panels show the results for a local NG with two values of /nl- 



orders in the Taylor expansion of the function in the exponential: 
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(D.6) 



Here we used the fact that ^'^'(A,) = 0{e) and ^("'^(A*) = 0{f^). The integrals in the second equality of this 
derivation can be computed analytically, which allows one to go to arbitrary accuracy with the saddle point 
technique. Notice that the results of these integrations are of higher order than the terms we retain. In the 
main text, where the integral contains also a polynomial 7^ (A) one can again compute the errors via similar 
Taylor expansions. These errors can be shown to be of order ©(e^), comparable to other terms which we 
ignore. 

One can also estimate the errors introduced by our opproximations by using the following toy model in 
which everything is computable: Take the 3-point cumulant ei to be different from zero and all higher order 
cumulants e„ for n > 2 to be zero^°. For such a model the integral is 



dA e*''-^"^^/2+(-'^)'''^i/6 
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1/2 



^/TT2e^ 



exp 



3il 



(D.7) 



In the r.h.s of this equation we have used the saddle point approximation but have made no expansion in ev. 
By comparing the numerical integration of the l.h.s. with the expression on the r.h.s. (panel (a) of Fig. 7), 
one can see that the errors introduced by the saddle point approximation are indeed of order e^ as indicated 
by (D.6). On the other hand, one can use the numerical integration of the left hand side of this equation and 



^"This toy model is inconsistent because if the third cumulant is different from zero, then all higher cumulants must 
also be different from zero. We use it here only to estimate how good the saddle point prescription is in approximating 
an integral, and compare it with errors induced by a perturbative expansion in ev. 
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compare it with the approximation (D.5) (panel (b) of Fig. 7), to see that the biggest error is of order i^v^ 
induced by the fact that we perforin a perturbative expansion in ev. Notice that here we considered only the 
leading order in tv and ignored unequal time correlators, while in the main text we present a result which is 
more precise (to next to leading order in tv) and complete (using the excursion set formalism rigorously). 
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